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We analyze a system of two-component fermions which interact via a Feshbach resonance in the 
presence of a three-dimensional lattice potential. By expressing a two-channel model of the reso- 
nance in the basis of Bloch states appropriate for the lattice, we derive an eigenvalue equation for 
the two-particle bound states which is nonlinear in the energy eigenvalue. Compact expressions 
for the interchannel matrix elements, numerical methods for the solution of the nonlinear eigen- 
value problem, and a renormalization procedure to remove ultraviolet divergences are presented. 
From the structure of the two-body solutions we identify the relevant degrees of freedom which 
describe the resonance behavior in the lowest Bloch band. These degrees of freedom, which we call 
dressed molecules, form an effective closed channel in a many-body model of the resonance, the 
Fermi resonance Hamiltonian (FRH). It is shown how the properties of the FRH can be determined 
numerically by solving a projected lattice two-channel model at the two-particle level. As opposed 
to single-channel lattice models such as the Hubbard model, the FRH is valid for general s-wave 
scattering length and resonance width. Hence, the FRH provides an accurate description of the 
BEC-BCS crossover for ultracold fermions on an optical lattice. 



I. INTRODUCTION 

In the absence of interactions, bosons cooled to zero 
temperature undergo a phase transition into a coherent 
quantum state of matter, the Bose-Einstein condensate 
(BEC), in which all particles reside in one single-particle 
state. Fermions, on the other hand, evolve smoothly from 
a classical gas into a degenerate Fermi gas where the 
N particles occupy the N lowest energy single-particle 
states. However, in the presence of arbitrarily weak in- 
teractions the Fermi gas has an instability towards a 
the Bardeen-Cooper-Schrieffer (BCS) superfluid in which 
fermions are bound to form bosonic Cooper pairs and 
these pairs then condense to form a BEC p] . In this case 
there is a transition from the normal gas to the BCS su- 
perfluid which occurs at temperatures much lower than 
the characteristic temperature which marks the onset of 
quantum degeneracy. As the attractive potential between 
fermions increases, the binding energy of a fermionic pair 
becomes much larger than the Fermi energy and these 
pairs become tightly bound bosonic molecules, a phe- 
nomenon known as the BEC-BCS crossover. The fact 
that this crossover is smooth implies that we can mean- 
ingfully speak of the crossover as Bose condensation of di- 
atomic molecules with coupling-dependent properties [2] . 
Thus, an understanding of these diatomic molecules, the 
two-body bound states, is key to understanding the be- 
havior of the system throughout the crossover. 

The BEC-BCS crossover has fascinated theorists for 
more than three decades [2HS], but it is only within the 
last few years that it has been amenable to experimen- 
tal study [BHTP]. One prevalent modern context of the 
BEC-BCS crossover is ultracold fermionic atoms in opti- 
cal traps interacting via a Feshbach resonance. In partic- 
ular, this work focuses on Feshbach interacting fermions 
in an optical lattice, a periodic potential made of in- 
terfering laser beams. The lattice introduces richness 
into the problem which is not encountered in the contin- 
uum or harmonic traps. For example, because of the re- 



duced translational symmetry, fermions can pair to form 
molecules with a center of mass quasimomentum differ- 
ing by a reciprocal lattice vector, leading to the possibil- 
ity of multiple bound states for a fixed s-wave scattering 
length. Additionally, the center of mass, relative, and 
internal degrees of freedom of an object in a lattice do 
not separate as in the continuum or a harmonic trap. 
This gives rise to bound state properties which depend 
both on the internal state and the center of mass mo- 
tion of the object. Lastly, new many-body phases can 
be induced by a lattice, such as the Ising deconfinement 
transition between an atomic superfluid and a molecu- 
lar superfluid recently predicted for Feshbach interacting 
bosons in lattices 

A prevailing theme in modern condensed matter 
physics is that the relevant degrees of freedom in a 
strongly interacting many-body system are often not 
the microscopic degrees of freedom. This realization is 
the underpinning of the Landau Fermi Liquid theory, in 
which quasiparticles with renormalized mass and interac- 
tions form a liquid in one-to-one correspondence with a 
nonintcracting Fermi gas. The renormalized parameters 
of such quasiparticle theories can often be determined 
from few-body physics by an appropriate dressing of the 
interaction via a Bethe-Salpeter equation or a renormal- 
ization group procedure [12J . In our situation the scatter- 
ing of two fermions at low energies can be well described 
by a few scattering channels in the continuum, but in 
the lattice each of these channels inherits a band index. 
The interchannel coupling is typically much larger than 
the band gap for a broad Feshbach resonance, and so 
all of the bands become strongly coupled [57]- This sug- 
gests that the bare channels which are relevant in the 
continuum are no longer an appropriate description for 
the strongly correlated physics in the lattice, as we must 
sum over a large number of bands for an accurate so- 
lution. Instead, we use the following facts to identify 
relevant dressed channels for the lattice system: (i) the 
interchannel coupling is large compared to the other en- 
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ergy scales of the problem; (ii) the interaction is short 
range compared to the other length scales of the prob- 
lem, in particular the lattice scale; and (iii) the density 
per unit cell is low. A table of scales of the lattice prob- 



lem is given in Sec. Ill 



Our approach to identifying the relevant degrees 
of freedom can be described in four steps, displayed 
schematically in Fig. [I] First, we identify the channels 
which are relevant to the scattering of fermions at low 
energies in the continuum. In this work we consider a 
two-channel model for a Feshbach resonance, but the ex- 
tension to multichannel models is straightforward. While 
we consider only two channels, corresponding to two in- 
dependent states in the continuum, these two channels 
become an infinite set of states in the lattice, indexed 
by band indices. The open channel is spanned by the 
scattering states of two particles, and so is indexed by 
two band indices. In one dimension (ID), the open chan- 
nel states in the lattice form scattering continua which 
are separated by band gaps as shown schematically in 
Fig. [T] The closed channel is a bound, point-like boson, 
and so is indexed by a single band index, with dispersion 
shown schematically in Fig. [I] The two-channel model 
parameters, the coupling constant g and the detuning is, 
which will be discussed in more detail in Sec. |IV| also 
inherit band indices. In this way all of the lattice chan- 
nels become strongly coupled when g is large. Second, 
we partition the Hilbert space into low energy and high 
energy parts, where low and high energy are defined rel- 
ative to the band gap. The low energy space consists of 
fermions in the lowest band of the open channel, and the 
high energy space consists of (a) all bands of the open 
channel which have at least one particle not in the lowest 
single-particle band and (b) all bands of the closed chan- 
nel. The high energy piece is only accessed via the Fesh- 
bach coupling, which is the largest two-body energy scale. 
Additionally, as the Feshbach coupling has an effective 
range much smaller than the lattice spacing, the high 
energy portion of Hilbert space is only accessed when 
two particles come in very close spatial proximity. At 
low lattice fillings, densities of more than two particles 
per unit cell are rare, and so we neglect interactions in- 
volving more than two particles. Accordingly, the third 
step in our effective model transformation is to numeri- 
cally solve the high energy piece for two particles. Our 
specific numerical method is discussed in Appendices [E] 
and[Cj The dressed eigenstates of the two-particle high- 
energy sector of Hilbert space form a more appropriate 
basis for describing scattering in a lattice at low energies 
than the bare bands. In the final step, we re-couple the 
low energy sector to the dressed high energy sector at 
the many-body level, resulting in a two-channel lattice 
model with renormalized parameters which we call the 
Fermi resonance Hamiltonian (FRH) [13] . 

There have been many other works, referenced and dis- 
cussed in detail in Section |Hj which have attempted to 
identify the relevant dressed degrees of freedom for lat- 
tice fermions. A key difference between our works and 
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FIG. 1: (Color online) Schematic of effective model trans- 
formation. We model the continuum scattering of fermions 
using a two-channel model for the Feshbach resonance. In the 
presence of the lattice, all channels inherit band indices. The 
Hilbert space is partitioned into low and high energy sectors 
relative to the band gap, and the high energy sector is diago- 
nalized at the two-body level. The two-body solutions of the 
high-energy piece are used as an effective closed channel and 
re-coupled to the low energy sector at the many-body level, 
resulting in the Fermi resonance Hamiltonian. 



others is that the solution of the high-energy piece, step 3 
in our scheme, is done using the full lattice potential and 
not a tight-binding or harmonic approximation. Approx- 
imations such as these which cause artificial separation 
of the center of mass and relative coordinates lead incor- 
rectly to both qualitative and quantitative differences in 
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the effective model when the interactions are strong. 

This article is organized as follows. Section [TT] re- 
views theoretical approaches to describing the BEC-BCS 
crossover in cold Fermi gases with a focus on the crossover 
in trapped geometries. Section III provides an overview 



schematic form 



of the physical scales of the problem. Section |IV| derives 
an equation for the bound states of two fermions in a 3D 
optical lattice interacting through a Feshbach resonance, 
discusses the symmetries of the solutions, gives explicit 
expressions for the matrix elements of the interchanncl 
Hamiltonian in the limit of an infinite lattice, and pro- 
vides details of the regularization of the theory. SectionfV] 



uses the results of Section IV to derive an effective many- 
body Hamiltonian for Feshbach interacting fermions at 
low energies and low density. We conclude in Section 
VI Some practical and numerical details concerning the 
two-particle solution are given in the appendices. 



II. OVERVIEW OF THEORY OF STRONGLY 
INTERACTING FERMIONS 



The simplest approach to the BEC-BCS crossover is to 
use the generalized BCS ansatz, which is valid only for 
weak interaction and high density, at the mean-field level 
for arbitrary interaction and find the chemical potential 
self-consistcntly by fixing the average number of parti- 
cles pQ. This approach was first used by Eagles in the 
context of superconductivity in low carrier concentration 
systems [3] and later by Leggett [5j and Nozieres and 
Schmitt-Rink [3] explicitly for the BEC-BCS crossover 
at zero temperature and finite temperature, respectively. 
The last two works demonstrated that the crossover is 
smooth, and that the transition temperature into the su- 
pcrfluid phases is smooth as a function of the coupling 
from the BCS to the BEC limit. Treatment at the mean- 
field level is quantitatively accurate for so-called narrow 
resonances where the width of the resonance is much 
smaller than the Fermi energy [14] , as this provides a di- 
mensionless small parameter in which to perturbatively 
expand. However, current experiments with ultracold 
atoms work in the opposite limit of resonance width much 
greater than the Fermi energy where the absence of any 
small parameter makes analysis difficult. 

Further work on the BEC-BCS crossover in the cold 
atomic gases context noted that the atom-molecule Fes- 
hbach coupling which generates the attractive interac- 
tion in alkali gases results in a condensate of molecules 
which is mutually coherent with the Cooper paired 
fermions [15] . This novel many-body state, called a res- 
onance superfluid, has a transition temperature near the 
divergence of the scattering length which is larger than 
that of the pure BEC or BCS limits. Models for reso- 
nance superfluids [15-20 include a pairing term of the 
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where ' (r) is a field operator for a fermion with spin 
a € {t, 1} and $>( b > (r) is a field operator for a bosonic 
molecule. Such a term converts a pair of fermionic atoms 
in opposite spin states into a molecule at the center of 
mass, and vice versa. In addition, to describe a Feshbach 
resonance one includes an energetic detuning hv between 
the open channel fermions and the bosonic closed channel 
molecules. To our knowledge the first time such terms 
were used in a cold atom context was in Ref. |21) . al- 
though such terms had been used phenomenologically in 
the study of high-T c superconductivity for many years 
prior [22-24]. The analogous high-T c model, known as 
the Cooperon model, is still a subject of current re- 
search [2"S] . 

The approximation that the pairing of two fermions 
to form a boson occurs only at the center of mass in- 
troduces an ultraviolet divergence in the theory. The 
proper renormalization of this divergence for the two- 
channel resonance Hamiltonian was first carried out in 
the absence of a lattice by Kokkelmans et al. [1T\. Here, 
two-channel means that only a single bosonic field cou- 
ples to the fermions. Kokkelmans et al. presented general 
relations between the bare and renormalized properties 
of the two-channel resonance model in terms of a momen- 
tum cutoff K* . Of particular interest is the renormaliza- 
tion of the detuning, in which the physical detuning v is 
related to the divergent bare detuning v as 



(2) 



with (?k the matrix element coupling two fermions with 
relative momentum k to a boson and et = h 2 k 2 /2m is 
the free fermion dispersion. For the contact interaction, 
.9k = 5, and the renormalization Eq. ^ becomes |17j 



hv = hug — 



mK*g 2 
2n 2 h 2 



(3) 



This renormalized detuning defines the zero-energy limit 
of the T-matrix [22 El], 
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Equation Q provides the correct behavior of the s- 
wave scattering length near a Feshbach resonance for any 
detuning, indicating that the renormalized two-channel 
model is an appropriate description of the resonance. In 
Sec. |IVD[ we discuss how to properly renormalize the 
lattice theory. 

It should be stressed that the low-energy limit of 
the two-channel resonance theory, Eq. Q, involves only 
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renormalized quantities. An alternative and more con- 
ventional formulation for finding the effective interaction 
between fermions near a resonance is to integrate out 
the bosonic field using a Hubbard-Stratonovich transfor- 
mation. However, the theory resulting from a Hubbard- 
Stratonovich transformation is expressed solely in terms 
of bare rather than renormalized quantities, and so the 
two approaches are inequivalent. Another means of inte- 
grating out the closed channel is replace the two-channel 
model with an effective single-channel model. In the 
single-channel model, a contact pseudopotential interac- 
tion between fermions is introduced, and is chosen to 
reproduce the proper scattering length. Such a pseu- 
dopotential causes difficulties near resonance where the 
scattering length diverges, as the T-matrix also diverges. 
In the two-channel model, the effective potential between 
fermions is the separable potential 
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and V denotes the Cauchy principal value [25J I2Z]- 
Hence, the two-channel model has residual momentum 
dependence in the effective potential even in the limit 
of a contact coupling g^ — » g. This momentum depen- 
dence makes the theory well-defined even at v = 0. The 
residual momentum dependence in the effective fermion 
interaction predicted by two-channel models is an impor- 
tant difference with single-channel models. 

In the presence of a lattice the natural choice of single- 
particle states are Bloch functions, and so an expansion 
of the field operators in terms of Bloch functions 

*(/)(r)=£ £ a nq ^(/ q ) CT (r) , (7) 

n qGBZ 

* (6) (r) = E E W^r) , (8) 

n q£BZ 

is appropriate. Here, (r) is a Bloch function diago- 
nalizing the fermionic single-particle Hamiltonian with n 
a band index and q a quasimomentum index, a nqCT cre- 
ates a particle in state (pnqa (r), </>n q (r) is a Bloch func- 
tion diagonalizing the single-boson Hamiltonian, 6 nq cre- 
ates a particle in state 0„ q (r), and BZ denotes the first 
Brillouin zone. Using these expansions in the Hamilto- 
nians describing resonance superfluidity leads to Fermi- 
Bose Hubbard Hamiltonians (FBHHs) of the generic form 
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tuning between the open and closed channels, E^q is 
the single-particle energy of a fermion in Bloch state 
<^£kj<t (r)j Eaq is the single-particle energy of a boson in 
state </>n q (r) , and "background terms" denotes spatially 
local interaction terms such as appear in the Hubbard 
model |28j arising from the background scattering 
length. 

Lattice FBHHs such as Eq. ^ generically require a 
large number of bands for an accurate solution. Because 
current methods used to solve models such as the FBHH 
scale poorly with the number of bands, many authors 
have attempted to derive effective models from a FBHH 
which are valid in some limiting case. For example, Carr 
and Holland [29] take a lowest band approximation of 
a FBHH, and show that in the limit of perturbative 
tunneling it becomes an XXZ model describing tunnel- 
ing of paired fermions. Zhou [30] begins with a similar 
lowest-band model, and demonstrates that Mott insu- 
lating states are unstable with respect to fermion-boson 
conversion at the mean-field level. Other authors do 
not take a single-band approximation, but rather group 
the excited fermion bands as well as the closed chan- 
nel molecules into a single dressed molecule field whose 
properties are fixed by two-body physics. We shall re- 
fer to such Hamiltonians as dressed effective Hubbard 
Hamiltonians (DEHHs). We will review several such ap- 
proaches below to contrast with the methodology taken 
in the remainder of this work. 

Dickersheid et. al. [31] consider a DEHH, and deter- 
mine the properties of the dressed molecule by consider- 
ing deep lattices and replacing the lattice with a single 
harmonic well. The harmonic trap approximation is the 
source of two systematic errors. First it artificially leads 
to separability of the center of mass motion from both 
the relative motion and internal structure. Second, it 
underestimates the extent of Wannier functions, often 
by an order of magnitude. We will show that qualitative 
properties of the tunneling, as well as its general order of 
magnitude, cannot be accounted for using this approach. 

The exact solution for two particles interacting via a 
Feshbach resonance in an anisotropic harmonic trap was 
determined by Diener and Ho [32]. In their calculations 
they stressed the importance of intra- as well as inter- 
band coupling terms, and properly renormalized the the- 
ory to remove divergences from using a point-like boson. 
Our results for the bound state energies reduce to theirs 
in the appropriate limit of deep lattices. However, their 
approach is still incorrect for determining the Hubbard 
parameters in a DEHH, as the harmonic oscillator ap- 
proximation suffers from the same systematic errors as 
Dickershied et al.. For example, the error in the tunneling 
is 70% for moderate lattice depths in use in experiments. 

Duan has derived DEHHs using both a projection op- 
erator formalism |33j and general symmetry consider- 
ations [Mj. In his approach, the dressed molecule is 
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determined by the solution of a two-atom, multi-band 
Schrodinger equation governing the atom-molecule con- 
version process on a single site. He then argues that 
when one of the dressed molecular energies is close to 
the scattering threshold of two particles in band n the 
only allowed on-site configurations are empty sites, sites 
with a single fermion in band n, and sites with a single 
dressed molecule. When this is the case, the Hamilto- 
nian can be projected into the reduced Hilbert space with 
only these on-site configurations. Because on-site pair- 
ing has already been taken into account in the definition 
of the dressed molecules, pairing in the effective many- 
body Hamiltonian occurs only between dressed molecules 
and fermion pairs which do not reside all on the same 
site. His work does not give a prescription for solving or 
renormalizing the on-site problem, but references Diener 
and Ho's work with the harmonic oscillator approxima- 
tion, which, as we stated, suffers from systematic errors 
in the calculation of the DEHH parameters. More recent 
work by Kestner and Duan [35] uses the numerical so- 
lution from a double-well potential to avoid some of the 
shortcomings of the harmonic oscillator approximation. 
While this work captures some of the physics of the lat- 
tice at the nearest-neighbor level, it does not capture the 
full quasimomentum dependence of the lattice Hamilto- 
nian. Other approaches [351 [37] have also taken partial 
account of anharmonic terms by series or perturbation 
expansions. 

Biichler was the first to give the exact solution for 
two fermions interacting through a zero-range Feshbach 
resonance in an optical lattice, properly accounting for 
the effects of higher bands and renormalization [35J . He 
then showed that when the interaction term U of the 
single-band Hubbard model was determined from the 
scattering properties of this exact two-body solution self- 
consistently the Hubbard model still failed to reproduce 
the correct physics even for moderate s-wave scatter- 
ing length. Specifically, the prediction of the two-body 
bound state energy from the self-consistent single-band 
Hubbard model deviates significantly from the exact re- 
sult for a s /a > 0.02 on the BEC side and for |a a |/et > 0.06 
on the BCS side [38j |39]. This has motivated our 
approach of using a two-channel lattice model instead 
of a single-channel lattice model such as the Hubbard 
model. As in our discussion of single-channel vs. two- 
channel models in free-space, the use of a two-channel 
model in the lattice allows us to capture the complete 
quasimomentum dependence of the scattering amplitude 
rather than only its low energy, low-momentum behav- 
ior. Buchlcr's discussion of the two-body solution focused 
on states with zero total quasimomentum, although the 
theory encompasses states with arbitrary total quasimo- 
mentum. In this work we show how extensions of his 
work may be used to derive a DEHH. 

Very recent work by von Stecher et. al. |40j focuses on 
a DEHH near a lattice resonance. Instead of solving an 
on-site problem they (i) project the two-body Hamilto- 
nian outside of the scattering continuum of two fermions 



in bands n and m which gives rise to the resonance, (ii) 
solve this projected Hamiltonian exactly, and then (iii) 
use the eigenstates of this projected Hamiltonian as a 
dressed closed channel. This approach is very similar in 
spirit to ours. However, von Stecher et al. assume low di- 
mensionality (quasi-lD) from the outset. Their approach 
breaks down when the energy associated with the Fesh- 
bach coupling becomes larger than the energy associated 
with the transverse confinement. In contrast, our model 
applies to arbitrary dimensionality. Our model treats 
the population in transverse excited states as being fixed 
by the two-body solution and thus part of the dressed 
molecule. Thus, any imposed condition of reduced di- 
mensionality, either quasi- ID or quasi-2D, can be con- 
trolled by the transverse tunneling and coupling rates of 
the dressed molecules. 

Finally, we note one other recent work by Titvinidze 
and coworkers [41] which accounts for the effects of higher- 
bands in a mean-field approximation assuming that the 
population of excited bands is small, and gives some pre- 
dictions of the resulting model using dynamical mean- 
field theory. 



III. ENERGY SCALES OF THE LATTICE 
FESHBACH PROBLEM 

To date, the most successful experimental realizations 
of two-component fermionic atoms exhibiting Feshbach 
resonant phenomena are 40 K [5j H2HS] and 6 Li [4"5"H4"5] . 
Typical operating temperatures are T ~ lOOnK in order 
to achieve quantum degeneracy at stable densities. In 
these systems, the two internal components of the atoms 
correspond to two sublevels Mp in a manifold of fixed 
total atomic spin F. Feshbach resonances [50] in these 
systems are due to hyperfinc couplings between scatter- 
ing states and a bound molecular state of the interatomic 
potential, and are tunable by magnetically shifting the 
position of the scattering threshold with respect to the 
bound state energy. Such resonances can be from milli- 
gauss to hundreds of gauss wide, depending on the species 
and state; in practice a well-controlled uniform magnetic 
field is sufficient to achieve a Feshbach resonance in ex- 
periments on many atomic species. The magnitude of the 
effective range R* of typical Feshbach resonances in these 
systems, which is related to the rs which appears in the 



low-energy scattering amplitude Eq. (131 in Sec. IV be- 
low as R* — —2ri,, is typically of order 20 angstroms or 
less [EE]. 

The lattice introduces a new length scale to the prob- 
lem, the lattice constant a. For optical lattices in the 
simple cubic arrangements that we consider in this pa- 
per, a = A/2, where A is the wavelength of the retro- 
reflected laser light forming the lattice. The energy scale 
associated with the lattice spacing is the recoil energy, 
E R = h 2 ir 2 /2ma 2 , which is of order 25kHz « 1/aK for 
40 K and 170kHz « 8/iK for 6 Li. The tunneling energy of 
a single particle in the lowest Bloch band along a princi- 
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pal axis of the lattice may be expressed as 
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where V is the lattice strength defined in Eq. (35). The 



coefficients in Eq. ( 10 ) represent the best fit to numerical 
data valid to < l%for V > 2Er. The tunneling gives 
a measure of the width of the lowest Bloch band. The 
gap between the lowest two bands in a simple cubic opti- 
cal lattice is nonzero for V > 2AEr, and monotonically 
increases with V/Er. Near V — 12Er, where the numer- 
ical examples in this work are performed, the band gap 
is approximately 5Er. The fact that the effective range 
is much smaller than the lattice spacing, rg/a -C 1 im- 
plies that the coupling parameter g which describes the 
interaction between channels in the two-channel model is 
large in the natural units of the lattice, see Eqs. ( 33p4 1 in 
Sec. |IV| below. The scales of the problem are summarized 
in Table U 



Scale 


Typical values 


Lattice spacing a (nm) 


550 


Temperature T (nK) 


-TOO 


Recoil energy h 2 iv 2 /2ma 2 (kHz) 


25 ( 40 K), 170 ( 6 Li) 


Lowest band tunneling t (Er) 


~0.01 


Band gap (En) 


~ 5. 


Effective range r_g (nm) 


£1 


Interchannel coupling g/Ena 3 ' 2 


>16 



TABLE I: Table of scales of the lattice Feshbach problem. 
The values given for the tunneling and band gap refer to V = 
12E R , see also Eq. JT61. 



IV. TWO-PARTICLE SOLUTION 

A. Derivation of the Nonlinear eigenvalue equation 

In this section we derive an equation for the bound 
states of two fermions interacting through a Feshbach 
resonance in the presence of a 3D optical lattice [T3] [35J • 
We begin by modeling the Feshbach resonance using a 
two channel model in the continuum. In this model, we 
partition our Hilbert space into open and energetically 
closed channels, with the asymptotic limit of the open 
channel potential corresponding to two free atoms. For- 
mally, this is achieved using a projector P into the open 
channel and its complement Q = 1 — P: 



(11) 
(12) 



where H is the complete two-particle Hamiltonian, \E' is 
the two-particle wavefunction, and, e.g., Hpp = PHP is 



(E - Hp P ) = HpqQ* , 
(E - H QQ ) = H QP P* , 



the piece of H taking the open channel Hilbert space to 
the open channel Hilbert space [52] • The closed channel 
potential Hqq is assumed to support a bound molecular 
state near the scattering threshold of the open channel 
potential. We take the open channel to be spanned by 
states of two fermions in different internal states with 
equal mass ra\ = mi — m and the closed channel to 
be spanned by molecular states with twice the fermionic 
mass and twice the fermionic polarizability [S3]. We de- 
note the two internal states of the open channel by the 
pseudospin notation {t, J,}. The two-channel model is 
defined by an interchannel coupling g which pairs the 
open channel fermions to a closed channel molecule at 
the center of mass and a detuning hv between the bound 
state of the closed channel and the scattering threshold of 
the open channel. The two-channel scattering amplitude 
is ED 



/(k) = - 
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(13) 



where [i is the reduced mass, k the incident momentum, 
and et the free particle energy. This scattering ampli- 
tude may be written in the asymptotic form of scatter- 
ing of particles with momentum small compared to the 
inverse effective range |54j by identifying the s-wave scat- 
tering length a s = —2/j,g 2 /8ir 2 ti i i' and the effective range 
tr = Trh 4 / fi 2 g 2 . Both a s and tr can be measured exper- 
imentally, for example by radio-frequency spectroscopy 
of Feshbach molecules [55] . 

Denoting the wave function of the two fermions in the 
open channel as ip (x, y) and the wave function of the 
closed channel molecules as (f> (z) , the two-channel model, 
Eqs. ( TT|T2 ), in position representation becomes 

[£-tf (x)-Ho(y)]V>(x,y) = 

g J dza(r)<f>(z)S(z-R) , (14) 
[E-hvQ-Hg 1 (z)]0(z) = 

g J dxdy a (r) $ (x, y) 5 (z - R) . (15) 

In this expression H (x) = — + Matt (x) is the sin- 

gle particle Hamiltonian for the open channel, H^ 1 (z) = 

~Tm^\ + 2Viatt ( z ) is the single particle Hamiltonian for 
the closed channel, a (r) is a regularization of the Fesh- 
bach coupling with cutoff A, and we have introduced the 
center of mass and relative coordinates R = an( j 
r = x— y. The subscript in Uq denotes that this is a bare 
detuning entering the microscopic theory which is related 
to the physically observable detuning v in the limit as the 
cutoff A — >• oo where a (r) — > S (r). We will discuss our 



particular choice of regularization in Sec. IV C 

The eigenfunctions of the single particle Hamiltonians 
are Bloch functions 



£ (x) 0<£ (x) 
?o M (z)<^(z) 



E%l<m (x) 



E. 



" q 



£(■) 



(16) 
(17) 
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which are described by a band index n = (n x ,n y ,n z ) 
and a quasimomentum q = (q x ,Qy,Qz) m the first Bril- 
louin zone (BZ). We establish the conventions that n and 
m are band indices for the open channel and their sums 
run from 1 to oo; s and t are band indices for the closed 
channel and their summations run from 1 to oo; q is a 
single-particle quasimomentum; K is the total quasimo- 
mentum; and sums over quasimomenta always run over 
the allowed values in the BZ. Because the interaction 
is invariant under translation of the coordinate system 
by a Bravais lattice vector, the total quasimomentum 
K = qi + q2 is conserved and we may diagonalize the 
Schrodinger equation in subspaces of fixed K. The open 
channel solution with total quasimomentum K may be 
parameterized as 

i>K (x,y)~EE ^ (q) <t>$ (x) <tf K - q W ■ 

* n in q 

(18) 

where iV 3 is the total number of unit cells in a 3D lat- 
tice with periodic boundary conditions. This wave func- 
tion generalizes the idea of single-particle Bloch states to 
multi-particle systems, as translating the system through 
a Bravais lattice vector R multiplies the wave function by 
a unimodular factor e — lK ' R . Using this form for the open 
channel wave function, parameterizing the closed chan- 
nel wave function as a sum over Bloch states in different 
bands as 



$ K (z) = £ T^W (z) , 



and inserting into Eqs. ( 14p5 ), we find 
[Ek - £* m (q)]vE» (q) = 

-j= v E ft "K (q) t s k , 

* S 

1 9 



(19) 



(20) 



(21) 



iV 3 ^ 



nm;q 



Here, we have made the definitions of the open channel 



energy E* m (q) = E, 
pling 



Kk (q) 

VWv 



(/) 

nq 



E. 



(/) 

m,K— q 



the interchannel cou- 



= J dxdy [<) (x) ^> K _ q (y )] * a (r) (R) . 

(22) 



and v, the volume of a unit cell. Note that we have affixed 
the subscript K to the energy eigenvalue to denote that it 
is the eigenvalue for fixed total quasimomentum. Taking 
the continuum limit allows us to make the replacement 



E 



iV 3 



VBZ Jbz 



(23) 



where vbz is the volume of the BZ. The equations ( 20]|2T ) 
become 



[Ek - E* m (q)]^ m (q) = 

* s 



(24) 



[E K -hv -EW]r?= (25) 

4:E/ — ^r(q)^ m (q), 
v v rr^ Jbz v bz 



Formally solving 
(E - H (x) - H (y) 



Eq. (24) by applying 



17] 



with rj a positive in- 



finitesimal and inserting into Eq. ( 25 ) gives 



[E K -v -E^]r? = 



(26) 



^E 



r/q 



E 



h*K (q) K£ * (q) 



BZ vbz t^Ejz- Ef m (q) + %n 



K 



The right-hand side of Eq. ( 26 ) diverges in the limit 



A — > oo, as is well known for two-channel theories involv- 
ing a point-like boson [T7]. We remove this divergence 
through renormalization, replacing the divergent bare de- 
tuning hvQ with the physical detuning hv by subtracting 
off the divergent part of the bare detuning [HI [26l [27] . 
Specifically, 



E^-hv- E^] T* = 9 - E xl (E K ) T* , (27) 

K (F , f dq \ - fegg (q) fefff* (q) K 
^bz ^bz ^E K - E£ m (q) + i?7 



K= y rfg v fegf (g) fe t n r (q) 
Jbz vbz -^ m (q) 



(28) 
(29) 



where the bars indicate that these quantities are com- 
puted in the absence of an optical lattice. The renormal- 
ization matrix xit ^ s a diagonal matrix, as will be shown 
explicitly in Sec. |IV D| The divergent parts of the two 
terms in Eq. ( 28 1 cancel each other and the limit A — > oo 
is finite. 

The renormalization prescription Eq. ( 28 1 is similar 



the one in free space, Eq. ([3|. There are two key differ- 
ences. The first is that the renormalization is not a scalar 
function as in free space, but rather a matrix indexed by 
closed channel band indices. The matrix element 
is the renormalization of the detuning of a closed chan- 
nel molecule in band s due to effective coupling through 
all bands of the open channel. The second difference 
is that, rather than imposing a high-momentum cutoff 
as in free space, we perform the regularization through 
the function a(r) The two procedures are the same if 
a(r) is chosen to be the Fourier transform of a Heaviside 
step function in the radial momentum. However, radially 
symmetric regularizations are inappropriate for the lat- 
tice renormalization due to a symmetry mismatch with 



the lattice. A better approach is to take a regularization 
function a(r) which is the Fourier transform of a func- 
tion with the symmetry of the reciprocal-space BZ. Our 
specific regularization and the manner in which the limit 
A — » oo is taken are discussed in Sec. IIVDI The com- 
putation of the overlaps h™{£ (q) is discussed in detail in 
Sec. |IV C| and efficient numerical methods for perform- 
ing the BZ integration in Eq. (281 and solving Eq. (27) 



are presented in Appendices [B] and [Cj respectively. 

The equation for the bound states, Eq. (27), is an 
eigenvalue problem which is nonlinear in the energy 
eigenvalue E^. The eigenvectors {T^-} satisfying this 
equation define the closed channel portion of the wave- 
function through Eq. ( 19 ). The open channel portion can 
be found using 



(q) 



9 \ ^ 



^"k (q) t£ 



V E K - E* m (q) 



IT] 



(30) 



together with Eq. (18). Here a is an index denoting in- 
dependent solutions of Eq. (27) for fixed g and v. Hence, 



the two-particle solutions of Eq. ( 27 ) may be written as 



*K Q (x,y) = ^f^T^(xHx-y) (31) 

C<- q (y)- 



T s K Q ^(q) 



where Afna is a normalizing factor and k(x — y) denotes 
a relative wavefunction for the closed channel which has 
characteristic volume v/A 3 . At the end of the computa- 
tion, we allow A — > oo, see Sec. |IVD| and so the purpose 
of k is to remind that the closed channel contributes only 
at the center of mass. The normalization coefficient is 



M£a = 1 - ( 
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E 



= ) 2 T 



K 



d X K (E£/E R ) 



dE 



T 



K 



K 



In the present work we specialize to the case of a sim- 
ple cubic optical lattice Matt (x) = yJ2 l= x, y ,z s ™ 2 (hi) 
with ki = it /a, a the lattice spacing. Thus, the volume 
of the unit cells in real and reciprocal space are v = a 3 
and i>bz = Stt 3 /v, respectively, and the natural unit of 
energy is the recoil energy E R — h 2 kf /2m. The parame- 
terization of the interchannel coupling and the scattering 
length in natural units arc 



(33) 
(34) 



9 


4 


j a 


E R a 3 / 2 


- n 3/2 


V r B 


a s _ 


2 

7T g 


Er 


a 


8 E 2 R a 3 


hv 



Additionally, the zero of energy is chosen to be the energy 
of two particles at zero quasimomentum in the lowest 
Bloch band, E° ± (0). 



B. Computation and symmetries of the 
single-particle basis 

In this section we discuss how to efficiently solve the 
single-particle Hamiltonians Eqs. ( 16p7 ) and discuss how 
the single-particle eigenfunctions transform under the 
point group symmetries of the lattice. The simple cubic 
lattice is separable in the sense that the lattice potential 
is a sum of ID lattice potentials 



Matt (X) = V 



=x,y,z 



sin' 



(hi) 



(35) 



Hence, the 3D Bloch functions </> nq (x) which form a com- 
plete eigenbasis of the single-particle Hamiltonian may be 
written as products of ID Bloch functions (j) nq (x) which 
satisfy 



P 

2m 



■ 2 ('> 

3m ( - 



4>nq (x) = E nq (j) nq (x) 



(36) 



Here, n denotes a one-dimensional band index and q a 
one-dimensional quasimomentum. We compute the ID 
Bloch functions numerically by writing 



<j) nq (x) = e lqx u nq (x) /VNa 



(37) 



where u nq (x) = u nq (x + a) are the eigenstates of the 
operator 



2m \ a J 



(38) 



Because of their periodicity, the functions u nq (x) may be 
expanded in a Fourier series as 



l 



i- 



u nq (x) = hm ]T c4e»/a . 

3 = ~l 



(39) 



(32) The real vector with elements {c J n } satisfies the eigen- 



value equation ^ ., H 



j. i 



— E r J 



where H^, 



-n,q J ^nq'-n,g 

(2.7 + q/hf /E R for j = f , -V/A for \j - j'\ = 1, and 
otherwise. The fact that the matrix H q is tridiagonal 
stems from the fact that the lattice potential, Eq. (35), 
contains only three Fourier components. Other lattices 
can be treated in an analogous manner. Convergence is 
achieved to machine precision even for high band num- 
bers ~ 20 by taking a finite Fourier cutoff / on the order 
of a few tens. 

In addition to translational symmetry, the single par- 
ticle Hamiltonian resulting from the simple cubic lattice 
is invariant under inversions of any set of Cartesian co- 
ordinates. Here, we discuss the transformations of the 
single-particle functions </> nq (x) in detail, as it provides 
a useful means to classify the interacting two-particle so- 
lutions, Eq. (31), see Sec. VB Additionally, proper care 



of the lattice symmetries is important to ensuring that 
the parameters in the effective many-body model trans- 
form sensibly under point group transformations. We 
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will denote the operator performing a spatial inversion 
as On, where 1Z is the set of Cartesian dimensions which 
are to be inverted. Additionally, an inverted coordinate 
will be denoted with a prime. That is, O-jz'k = x' where 
x'j = —xj for all j € TZ and x'j = Xj otherwise. When 
the Hamiltonian is time-reversal invariant the presence 
of an inversion symmetry generated by 6-r, implies that 
the Bloch functions of the simple cubic lattice transform 



0nq (x') 



II( 

j en 



0nq' (x) 



(40) 



At the extremal points of the BZ, q = and q = n/a, 
translations and inversions commute and so the associ- 
ated Bloch functions have well-defined parity with re- 
spect to inversion. At all other points in the BZ, trans- 
lation and inversion do not commute, but the inversion 
operator connects degenerate states at different points 
in the BZ up to a phase. From Eq. (40), we see that 
u nq (— x) = (—l) n+1 u n _ q (x). For the extremal points 
in the BZ this implies that the vectors c nq must have 
definite parity about j = and j = 1, respectively, see 
Eq. ( 39 1 . When computing the vectors c nq numerically 



the even and odd parity states become quasi-degenerate 
for high bands where the band gap becomes smaller than 
machine precision and the numerical eigenvectors will not 
automatically have definite parity. However, it is always 
possible to sort the eigenstates according to their parity 
to maintain the relation u nq (— x) — (— u nt - q (x). 
There is still an additional global phase ambiguity be- 
tween Bloch functions with different values of the quasi- 
momentum q. We fix this phase ambiguity by requiring 
that the Bloch functions are smoothly varying functions 
of q. This phase prescription gives rise to maximally 
localized Wannier functions [5B] , 



The inversion relationship Eq. (40) implies that the 



matrix elements of the inter-channel coupling transform 
as 



K£ (q') = n ( 

jen 



-1 



rij+mj+Sj+l j nm 



and so 



xft (jsk) = n (- r 

jen 



K£ (q) 



(41) 



(42) 



At extremal points of the total quasimomentum BZ, 
K = and — fc; (1, 1, 1), the inversion properties of x K 
imply that only molecular bands which transform iden- 
tically under complete inversions mix. Numerically, we 
find that the mixing between bands even for arbitrary 
K is generally small for energies near the lowest open 
channel band except at isolated points where eigenval- 
ues cross. Inversion leaves the eigenvalues of x K (Ek) 
invariant, but the eigenvectors transform as 



jen 



(43) 



where P a denotes a unimodular phase. Because the in- 
teraction is invariant under all elements of the inver- 
sion group we can classify the eigenvectors accord- 
ing to their inversion properties [13] . Accordingly, we 
choose P a — rijeK (— l) Sj+ , and require that be a 
smooth function of K otherwise. This fixing of arbitrary 
phases implies that the interacting two-particle solutions, 
Eq. (31), are also smooth functions of K. Importantly, 



this convention also ensures that the Wannier functions 
have their expected behavior with respect to the symme- 
tries of the lattice, see Sec. 



VB 



C. Matrix elements of the inter-channel coupling 

We now focus on the computation of the matrix ele- 
ments of the inter-channel coupling, Eq. ( 22 ) . We take 
as our regularization 



a(r) 



dk 

v(A) (27r) 3 



-ikr 



(44) 



where v (A) = .vbzA 3 and A is a positive integer. This 
regularization is similar to the momentum-shell regular- 
ization with cutoff K* in free space, see Eq. pi). However, 
rather than choosing a spherically symmetric cutoff, we 
choose the cutoff function to have the same symmetry as 



the BZ. The integral Eq. (44 1 may be explicitly evaluated 
to yield 



a r 



n 



sin (Atttv) 



(45) 



which is separable along principal axes and an even func- 
tion. For the simple cubic lattice, this implies that the 3D 
overlaps Eq. (22) are products of ID overlaps. Objects 



in ID will be distinguished from their 3D counterparts 
by not having bold symbols as arguments. We note that 
each overlap as written is dimensionless. 
The ID overlaps take the form 

j-N j-N 

Kk (q) / d Xl / dx 2 ( Xl ) £] K _ q (x 2 ) 

Jo Jo L 

x ^^i+^y {xi _ X2h (46) 

The highly oscillatory nature of both the Bloch functions 
and the regularization makes the numerical evaluation of 
this integral by quadrature computationally expensive. 
Rather, we seek to find an analytic expression for this 
integral in terms of the Fourier components of the Bloch 
functions, see Eq. (39). 



We begin by changing integration variables to £ = 
(xi + X2) /2, rj = (xi — X2) /2 with Jacobian 2. The re- 
gion of integration is bounded by the parametric (77, £) 
curves as (x,x), x £ [0, ^1 j ( x t~ x ), x £ — t]> 
(x,N-x), x e [f,0]; and (x,N + x), [-f,0]. The 
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integral becomes 

h fM= E p d ^v)f d^Uo 
2ViV P ={-i,i} Jo 



i" i 



(47) 



Expanding the Bloch functions as in Eq. ( 37 1 , we have 
(q) 



E P 



pN/2 



dr\e 



-ir)(2q-K-2irT) 



a (2t?) 



(i 



p={-i»i> 

pN—pr} 

< I dHe-^u sK (0 «*, (£ + 77) « - r?) , 

Ipri 

(48) 



where T is an integer such that 2irT shifts K — q to lie 
in the BZ. If we now expand the functions u in terms of 



their Fourier components according to Eq. ( 39 1 , we find 

2 = ££ E E K^.K-^ 



j,j'J"=-l p={-l,l} 



pN/2 



x / d7 7 e- i "[^- K ) +2 ' r ^- i '- r )]a(2r7) 





pN—pr) 

/ ^-^[r+j+j'-;"] 
Jpij 



(49) 



Changing integration variables to 77 = prj and using the 
fact that a (77) is an even function yields 



j,j'J"=—l 

x / d?7 cos (27Tz?7) a (2t?) / d£e i27r£t . 



where 

*=^^ + (j-j -r) 
^i"-i-/-r. 

We now focus on the £ integral: 



(50) 

(51) 
(52) 



JV-7) 



£t= I ^(e^-^-e 2 "**) f^O 
I N — 2-q otherwise 



Noting that t is always an integer, we have 

N ~\e^=N8 tfi -^^. (53) 

Tit 

Hence, the overlaps may be written 

l 

^( 9 )= lim2 E < g 4l, K _ ? ci+/ +r 7 1 (iV,A ) z) 



j,j'=-l 



3,j',j"=-l 



d K I 2 (N,A,z,t), (54) 



where 



7 1 (JV,A,z)= f 7 dr/ Sm(27rA7 ^ cos(27T77z) , (55) 
r ,»r . \ f^ 2 , sin(27rA77) r „ , sin27TTrf 

I 2 (N,A,Z,t)= dn ^ ^C0S27T77Z -. 

Jo T"! nt 



(56) 



In these integrals we have substituted the actual form of 



the regularization, Eq. ( 45 ) . The first integral is 



h(N,A,z) = — [Si (Mr (A + z)) + Si(A%(A-2))] 

Z7T 



(57) 



where the sine integral Si has the definition 



Si (y) = dx 



SIM 



(58) 



For the second integral, we first note that 



A (cos (A7rA) cos (Nttz) - 1) , . 
Kmh(N,A,z,t)= { \ 2{ ,2_ A 2 ) } ( 59 ) 



As we anticipate taking the limit A -> 00, we may choose 
that A is even without loss of generality and so 



,. -r ,*r . \ A (cos (A7rz) - 1) 

lim h{N, A, z, t) — y \ > ' , 60 

t->0 TT Z [z* — A Z ) 



which vanishes as z — > A. For t =/= 0, we have 

4n 2 tI 2 (N, A, z,i) = C (N, t - z - A) + C (N, t + z - A) 
- C ( N. t - z + A) - C (A, t + z + A) , 

(61) 

where 



C(N,x) = Ci (Attic) -logx, 



(62) 



Ci(a:) = 7 + loga;+ / dt— — -. (63) 



Here, 7 w 0.5772 is the Euler-Mascheroni constant and 
Ci (x) is the cosine integral. The function C (A, x) is even 
in x, and so we may take the absolute values of the x 
arguments without loss of generality. Furthermore, while 
the both the cosine integral and the logarithm diverge as 
x -> 0, lim^o C (A, x) = 7 + log Att. 

In summary, we have 
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Kk (g) — ~ c n,g < ^i,K-5 C s 



+j +r 



Si (TVvr (A + z)) + Si (A/tt (A - z)) + 



2A (cos (A/ttz) - 1) 



(64) 



Ntt (A 2 -z 2 ) 

E < S ( K -,4i [W I* - z - A|) + C(A/, \ t + z - A|) - C(AT, |* - « + A|) - C(iV, |t + z + A|)] , 



33' 3" 



where the prime on the latter summation indicates that 
terms where t — are excluded. Let us now con- 
sider taking the limit of an infinite number of unit cells, 
N — > oo. For the primed summation, we have that 
lim. E _j. 00 Ci (x) = and the logarithms are independent 
of N, so this entire summation vanishes as N — > oo, 
x 7^ 0. The only other possibility is that one of the 
terms in the absolute values is identically zero, but in 
that case lim__>.oC (N,x) = 7 + log Ntt. All such terms 
in the primed summation vanish as N — > 00 because of 
the factor of N in the denominator. Hence, the primed 
summation has no contribution in the limit of an infinite 
lattice. 



D. 



For the first summation in Eq. (64), we note that the 



last term in brackets vanishes as N — > 00, as z is real. 
When z = A this quantity vanishes identically for any 
N, as A is defined to be an integer. Hence, the only 
terms that remain are the sine integrals. We note that 
lmXj._5.00 Si (x) = 7r/2, but also that Si is an odd function. 
Thus, 



lim — 

N— >oo TT 



Si (Ntt (A + z)) +Si (Ntt (A - z)) 



-A < z < A 

1*1= A , 
otherwise 

(65) 
and so 

^lim h?$ (q) = E C n ^ K _ q ^/ +r (66) 
-2q-K + 2n (j-j'-vy 



33' 

x rect 



2ttA 



where rect (x) is the rectangle function. In the limit of 
an infinite lattice, the effect of the regularization is to 
cut off the contribution of the overlap in the shifted rel- 
ative Fourier space defined by z, see Eq. (51 1, with a 



rectangle function of width 2A. For finite N, the over- 
laps Eq. ( |64[ ) display a pronounced Gibbs phenomenon 
which make integration difficult near z = A. By tak- 
ing the N — > 00 limit the function becomes much better 
behaved. Additionally, taking N — > 00 allows us to con- 
sistently make the claim that the quasimomentum q is a 
continuous variable and so replace all sums by integrals 
according to Eq. (23) in Sec. IV 



Regularization of the theory and evaluation of 
the renormalization at K = 



In this section we use the results of Sec. IIV CI to find 
an expression for the renormalization matrix, Eq. ( 29 1 , 



at zero total quasimomentum. In addition to the theoret- 
ical appeal of being able to compute the renormalization 
analytically, this procedure also sheds light on the proper 
regularization procedure. Let us first define a shell sum- 
mation over bands with shell parameter S, 2~2 nm .g, as the 
summation over all n and m whose components are less 
than or equal to S with at least one of the components 
being S. To illustrate this notation, the shell summation 
for S = 1 consists of n = (1, 1, 1), m = (1, 1, 1): 



nm;l n x — 1 n y — 1 n z — 1 ra x — 1 m y — 1 m z — 1 

The shell summation for S = 2 is 

2 2 2 2 2 2 

E^EEEE E E 

rim: 2 n x — l n y — l n z —l in x — l m y — 1 m z — 1 
1111 1 1 

-EEEE E E • 

n x — 1 rin,— 1 n z — 1 m x — 1 m u — 1 m z — 1 



(67) 



(68) 



The subtracted summation ensures that at least of one 
the elements of n or m is 2, as required by our defini- 
tion. Shell summation provides a way to systematically 
add contributions from higher bands to a summation, 
as the summation over all bands with components less 
than or equal to £ is equal to the summation over all 
bands with components less than or equal to (t — 1) plus 
Snm f By induction, the summation over all bands with 

components less than or equal to S is __nm f ■ 

To properly renormalize the theory and to include the 
effects of all higher bands one must take both the shell of 
open channel bands S —¥ oo and the regularization cutoff 
A — s- oo. However, as shown in Ref. [3H], the order of 
these limits is important, and the proper limiting proce- 
dure is first to take the shell parameter S — ¥ oo with A 
held fixed and then let A — > oo. The limit of infinite cut- 
off A is taken using the asymptotic scaling relation 



[x5(2_0] (A)= Pst /A + x ? t (E K ), 



(69) 



where p s t is the slope defining the scaling with A and 
Xst(-^K) is the value as A — > oo. 
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Let us now consider computing the renormalization 
matrix, Eq. (29 1, at K = 0. Because the renormaliza- 
tion is computed in the absence of an optical lattice, the 
Fourier expansions of the single-particle states, Eq. (391, 



contain only a single plane wave, and the open channel 
energy becomes E® m (q) = 2Er<\ ■ q/7r 2 . The fact that 
all Fourier expansions contain a single plane wave im- 
plies that the overlaps h^(q) with s and K fixed are 
nonzero only for a specific pair of open channel band 
indices (n, m). This in turn implies that the renormal- 
ization matrix % K , defined in Eq. (29 1, is diagonal. The 



non-vanishing Fourier component for the lowest closed 
channel band with s = (1, 1, 1) is cjjn = Sj^. Further- 
more, the shift r introduced in Eq. (48 ) is when K = 0, 
and so 



hit 



(q)= II [e(g, + A)-e(g„-A)] , (70) 

u={x,y,z} 



with Q (x) the Heaviside function. Thus, provided that 
A > 1, we have that the 5=1 shell contribution to the 
s = 1, t = 1 element of the renormalization is 



relevant integrals are 



h =- 



1 



16E T . 



-1 ,2 

dx + dx 

2 Jl 
1 



dy 



x / dz-„ 

_ 1 x^ + y z + z z 



h = 
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Specifically, we have that Xix s=2 = ^-Tt + 3/2 + ^3- How- 
ever, we note that this is 
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(73) 



provided that A > 2. Following this construction for 
higher shells, we have the result that 



v° 

Xll;S=l 
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WEr 



dx 



dy 



dz- 



(71) 



This integral cannot be evaluated analytically, but its 
numerical value is w (-0.959265 ± 2 x \Q- e )/E R . The 
evaluation of this integral and the error estimation were 
performed using the Genz-Malik algorithm, discussed in 
Appendix [B] 



Let us now consider contributions from the second 
shell, 5 = 2. From the constraint that only j + f = 



terms contribute, see Eq. (66), we have that the only 
nonzero h™j) (q) overlaps along each dimension have n 
and m either both in the lowest band or both in the first 
excited band. Because of the shell summation constraint 
that at least one of the band indices is greater than 1, for 
each nonzero contribution we will have at least one pair 
of open channel band indices (n, m) in the first excited 
band. Now, integration over the energy of the excited 
band can be done in the extended zone scheme by inte- 
grating from q — —2 to — 1 and then from 1 to 2. The 



min(5,A) x°i-s=i < 



(74) 



In Sec. |IV C| we showed that the regularization param- 
eter A may be interpreted as a cutoff in Fourier compo- 
nents, see Eq. (66). The present analysis demonstrates 



that A can also be interpreted, albeit more loosely, as a 
cutoff in the contributions from higher bands. As the lat- 
tice potential vanishes, higher bands become equivalent 
to plane waves of high energy. Hence, the given lattice 
renormalization procedure produces the same results as 
the free-space renormalization procedure, Eq. ([3|, where 
the cutoff K* ~ An/a. For the first term in the defi- 
nition of \ K {Em) which includes the optical lattice, see 
Eq. (28), A is not a strict cutoff for the contributions 



of higher bands because the expansions Eq. ( 39 ) do not 
contain only a single Fourier component. This is why it 
is important to first take the limit 5 — » oo with A fixed 
and then take A — > oo. However, the contributions to the 
integral from shells 5 > A quickly approach their free- 
space values, causing the shell summation of x K (En) to 
converge, see Fig. [2] 



V. THE FERMI RESONANCE HAMILTONIAN 

We now proceed to step 2 of the scheme outlined in 
Fig. [l] the partitioning into low and high energy spaces. 
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FIG. 2: (Color online) Scaling of the effective lowest-band 
closed channel T-matrix x¥i with open channel band cutoff 
S and regularization cutoff A. (a) The T-matrix element 
Xii(-Ek) computed at bound state energy £"k = — IEr for a 
variety of open channel band cutoffs 5* and regularization cut- 
offs A. The open channel band cutoff S is increasing towards 
the bottom of the figure. The crosses are data and the solid 
lines are fits to the S — ¥ oo scaling form Eq. \69\ , with the 
solid red line being a fit through the first three points and the 
dashed blue line being a fit through all five points. The pre- 
dictions from these two fits differ only by a few percent, (b) 
The T-matrix element Xii computed at A = 6 (red pluses), 7 
(green crosses), 8 (blue asterisks), and 9 (purple boxes) for a 
variety of open channel band cutoffs S. Near S = A a marked 
change in behavior occurs followed by rapid convergence to 
the S — > oo limiting value. 



Typical operating temperatures of ultracold atomic gas 
experiments are below the band gap of the confining lat- 
tice potential, see Sec. |III| For moderate lattice strengths, 
this implies that fermions remain in the lowest open chan- 
nel band when they are separated by a distance larger 
than the effective range tb <C a to minimize their en- 
ergy, and this defines the low-energy, long-wavelength 
sector of Hilbert space. Only when fermions come to- 
gether at very close distance comparable to the effective 
range are they able to populate excited bands of the open 
channel or transfer population to closed channel bands. 
At low lattice fillings, the particular population of the 
high energy sector of Hilbert space is set by the two-body 



physics. We find the relevant states from the high-energy 
sector of Hilbert space by solving the two-body problem 
projected into the high-energy subspace. We call the 
eigenstates of the high-energy sector dressed molecules. 
The effective many-body description is then a resonance 
model between unpaired fermions in the lowest band and 
dressed molecules. By the construction of the dressed 
molecules, this model reproduces the correct scattering 
and bound state properties at the two-particle level, and 
hence is expected to provide an accurate many-body de- 
scription at low densities. 

We refer to the resulting model as an effective model in 
the spirit of effective models such as the t — J model [57] 
in that we have removed high-energy degrees of freedom 
and kept only the most relevant couplings for the low 
energy physics. While excited bands are still present 
in the theory in the form of the dressed molecules they 
are restricted to occur only in specific combinations, just 
as quarks are required to be bound at low energies. In 
solid state systems, often little is known about the mi- 
croscopic Hamiltonian, and so the construction of an 
effective Hamiltonian may not follow a well-controlled 
prescription. In contrast, the derivation of our effective 
model follows in a step-by-step manner from the micro- 
scopic physics. Our partitioning into low and high energy 
sectors is accomplished using projectors L into the lowest 
open channel band and D = 1 — L into all excited open 
channel bands as well as the closed channel molecules. 
A similar approach was taken in Ref. [40] . An analysis 



parallel to Sec. IV gives a nonlinear eigenequation for the 
closed channel components of the high-energy sector as 



[Ex-hv-Eg] T K 



9 ,~,K 



t 



X* (23k) Tf , (75) 



Xf t (E K ) 



E 



dg fex(qR n K*(q) 

vbz E K - E^ m (q) + ir] 



Xft 



(76) 



where the prime on the summation denotes that the sum 
is taken over all (n, m) ^ (1, 1). That is, the sum does 
not include the lowest open channel scattering band. The 
renormalization matrix % K appearing in Eq. ( 76 ) is the 



same as in Eq. (28), and so the detuning and scattering 



length appearing parametrically in this projected model 
are the renormalized detuning and scattering length of 
the full two-body problem using all of the open channel 
bands. In this way the projected system reproduces the 
correct, renormalized physics at the two-body level when 
recombined with the lowest scattering band of the open 
channel. 



A. Properties of the projected model 

The bound state energies of the non-projected two- 
particle problem were displayed as a band structure in 
Ref. |13j where the Fermi resonance Hamiltonian was 
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FIG. 3: (Color online) Comparison of the zero- 

quasimomentum bound states of the two-body solution and the 
projected two-body solution. All results are for a 3D cubic op- 
tical lattice of strength V = 12Er, and the grey stripes denote 
scattering bands of the open channel. The red crosses repre- 
sent the two-particle bound states with completely even parity 
under inversion, and the green asterisks have parity (1, 1, —1) 
and cyclic permutations. The solid red lines are the bound 
states of the projected model with completely even parity and 
the green dashed lines are the bound states of the projected 
model with parity (1,1,-1) and cyclic permutations. The 
odd parity solutions of the projected and full models very 
nearly coincide, but the resonance behavior in the completely 
even parity solution is drastically altered by the projection. 



first identified. In this section we examine the bound 
states of the projected two-particle problem for compar- 
ison. For clarity of exposition, we will study the prob- 
lem only at zero total quasimomentum. The results of 
the comparison are shown in Fig. [3j which shows the 
lowest bound state energies of the two-particle problem 
both with and without projection as a function of the s- 
wave scattering length in the limit of an infinitely broad 
resonance, see Appendix [C] The red crosses and solid 
red lines represent solutions with completely even par- 
ity under inversion p = (1, 1, 1) of the non-projected 
and projected models, respectively. The green crosses 
and dashed green lines represent solutions with parity 
p = (1, 1, —1) of the non-projected and projected models, 
respectively. The results for the odd parity bound states 
are almost identical between the non-projected and pro- 
jected models, which is a statement that this bound state 
does not couple strongly to the lowest open channel scat- 
tering band. In contrast, the lowest energy bound state 
with completely even parity displays drastic differences 
between the projected and non-projected models. The 
resonance position of the projected model, which is the 
energy with respect to the bottom of the lowest open 
channel scattering continuum where the s-wave scatter- 
ing length diverges, is shifted close to zero, indicating 
that this bound state causes a two-body resonance when 
energetically close to the threshold of the lowest band of 



the open channel. This is the sense in which the dressed 
molecules form an effective closed channel for the two- 
body resonance near the lowest scattering band of the 
open channel. 

As discussed in Ref. [T3] , the asymptotic scaling of the 



matrix \ with the regularization cutoff Eq. ( 69 ) holds 



for the projected matrix x (-E-k)i a s well. This follows 
from Eq. (66), which demonstrates that as A becomes 



larger than the support of the Fourier expansions of the 
lowest open channel band Bloch functions the overlaps 
^Ik (*?) are n0 longer functions of A. Hence, the scaling 
with A is dominated by the contributions of higher bands. 
The same argument extends to a projected matrix % K in 
which any finite number of open channel bands have been 
projected out. These results reinforce the interpretation 



put forth in Sec. IV D that the effect of the regularization 
is a cutoff in the contributions from higher bands, and 
so the results for large A should be insensitive to the 
behavior of low-energy bands. 

We conclude this section by noting that the method 
of projection could be applied in principle to some 
other scattering continuum (n, m), in which case the 
bound state which couples most strongly to that 
continuum would have its resonance shifted such 
that the scattering length diverges at approximately 
min {E* m (q) , q G BZ, K € BZ} . This allows for the 
derivation of effective Hamiltonians for resonances in ex- 
cited bands. Additionally, one can consider the first n 
scattering bands of the open channel to be the low energy 
sector of Hilbert space and project these out of the two- 
particle theory. The dressed molecules of the projected 
high-energy sector of Hilbert space will define an effective 
closed channel which reproduces the scattering properties 
of the lowest n open channel bands properly. In this way, 
the construction can be extended to higher relative scat- 
tering energy where scattering states in higher bands can 
play a role. In the remainder of this document we focus 
only on the case of resonance with the first scattering 
continuum. 



B. Wannier functions and Kohn construction 

We now move on to item 3 of the scheme outlined in 
Fig. [l] the solution of the high-energy sector of Hilbert 
space. The details of the numerical solution are discussed 
in Appendices |B||C| In this section, we focus only on the 
structure of the solutions. The two-particle solutions of 



the projected nonlinear eigenvalue problem Eq. ( 75 ) may 
be written as 



t K „(x, y) = -L[£ Tf a ^( X )«(x - y) 



9 

nms;q 



(77) 



£! K - q (y) 



where the prime on the summation indicates the exclu- 
sion of the lowest two-particle band of the open channel, 
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compare Eq. (31), and the reduced Feshbach coupling g 



was denned in Eq. (33). Here, the normalization factor 
is 



Kc 



k dx*(E£/E_ 



R) 



8Ei 



T 



K 



K 



(78) 



We define two-particle Wannier functions centered 
around some lattice site i as the quasimomentum Fourier 
transforms of the two-particle solutions, 



W ia (x, y) = 




-iKRi 



*Ka(x,y) 



(79) 



As is known from the seminal study of Kohn in ID [55] , 
the choice of phases for the Bloch functions affects the lo- 
calization properties of the Wannier functions around the 
site i. In ID, choosing the Bloch functions to be smooth 
functions of the quasimomentum leads to exponential lo- 
calization of the Wannier functions. In the present case, 
the phase s are chosen according to the construction of 
Sec. IV B in which the eigenvectors {T^} are smooth 



functions of the total quasimomentum except possibly at 
high-symmetry points of the BZ and phases under inver- 
sion are specified by the parity of the eigenstate. Ac- 
counting for transformation under inversion is important 
to ensure that the parameters in the resulting many-body 
model have the same symmetries as the few-body physics. 
The parity of the eigenstate is unambiguously defined at 
extremal points of the BZ in which all components of the 
total quasimomentum are either the zone center K = 
or the zone boundary K — —ir/a, and then extended into 
the general BZ by the smoothness requirements on the 
eigenvectors T^,. We will refer to this choice of phases as 
the Kohn construction. More localized functions may be 
possible by replacing the quasimomentum Fourier trans- 
form in Eq. ( 79 1 with a more general unitary transfor- 



mation which is consistent with inversion symmetry; we 
do not consider this possibility further here. 



C. Statement of the Hamiltonian and evaluation of 
parameters 



With the dressed molecule Wannier functions Eq. ( 79 1 , 



we can state the Fermi resonance Hamiltonian (FRH) 
as [13] 

-EE*sA4.+ ■£*«£*£ 

a£A4 i,j a£A4 i 

+ £ £fl£-k,k-i [d\, a a 3 ta ki + h.c 



(80) 



a£A4 ijk 



Here, is a fermionic operator which destroys a parti- 
cle with spin a in the Wannier state Wi a (r) constructed 
from the lowest band Bloch functions of the open chan- 
nel, n^J — a\ a ai a is the number operator for particles 



with spin a in the lowest band of the open channel, di a 
is a bosonic operator which destroys a dressed molecule 
in Wannier state W ict (x,y), and nfj is the number op- 
erator for dressed molecules in state a on site i. The 
summations over a £ M. indicate that the particular set 
of dressed molecules M. which are relevant to describing 
the resonance properly depend on the open channel band 
being projected, symmetry considerations, and the ener- 
getic level of approximation used to truncate the terms 
appearing in Eq. (80). For the case we focus on here 



in which the lowest open channel band is projected out, 
the lowest energy dressed molecule which has completely 
even parity under inversion is the most relevant one, as 
all others either have vanishing on-site effective pairing 
by symmetry or are far detuned from resonance [13j . 

We now turn to the last part of the program outlined 
in Fig. [T] in which the low-energy and high-energy sectors 
of Hilbert space are re-coupled at the many-body level. 
This re-coupling is performed by taking the expectation 
of the Hamiltonian within the basis of the open channel 
fermions in the lowest band and the dressed molecules, 
and provides us with the Hubbard parameters appearing 
in the FRH, Eq. J80|. The first term in Eq. p&) repre- 



sents tunneling of fermions with spin a in the lowest open 
channel band between neighboring sites i and j. The cal- 
culation of the tunneling using the Wannier functions for 
the lowest band fermions is well known [58 . The next 
term is the energy offset of a fermion in the lowest band 
with respect to the lowest two-particle scattering contin- 
uum, where Eq — J dq£'i,q/i>BZ- The remaining three 
terms represent processes involving dressed molecules. 
The first represents tunneling of the dressed molecular 
center of mass between two lattice sites i and j. The two 
sites i and j need not be nearest neighbors due to the non- 
trivial dispersion relation of the dressed molecules. The 
next term represents the effective detuning of a dressed 
molecular Wannier function away from resonance with 
the lowest two-particle scattering band of the open chan- 
nel. The final term is the effective pairing of two open 
channel fermions into a dressed molecule at the center of 
mass. 

The dressed molecular tunneling Hubbard parameter 
is evaluated as 

t?j = - / dxdyWt, (x,y)HW ja (x, y ) , 



dK 

VBZ 



iK-(Ri-Htf) F a 



(81) 



where H is the two-particle Hamiltonian. The fact that 
the tunneling does not change the internal state a of the 
dressed molecules is a consequence of the fact that the 
dressed molecules diagonalize the high-energy sector of 
the Hamiltonian. The remarkable fact that the tunnel- 
ing of these very strongly correlated objects is a simple 
quasimomentum Fourier transform of their dispersion is 
due to the invariance of the entire system under transla- 
tion. 
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The dressed molecular tunneling has several novel fea- 
tures compared to single-particle tunneling. The interac- 
tion destroys the separability of the band structure along 
principal axes and so tunneling is allowed along directions 
which are not the principal axes of the lattice. This is in 
stark contrast to single-particle tunneling, which always 
occurs along the principal axes. Additionally, for traps 
where the center of mass and relative motion separate, 
diagonal tunneling will always vanish. Thus, diagonal 
tunneling is a unique feature of the FRH due to non- 
separability in the dressed closed channel. While we have 
never found the diagonal tunneling to be greater than or 
equal to the nearest neighbor tunneling, it is often as 
large or larger than the second neighbor tunneling and 
in some cases can be of the same order of magnitude as 
the nearest neighbor tunneling. Second, the dispersion 
relation of the dressed molecules is highly nontrivial due 
to the fact that dressed molecules consist of two strongly 
interacting particles, and so the sign of the tunneling 
depends on the parity of the eigenstate and the direc- 
tion and magnitude of the tunneling in a nontrivial way. 
Quantitative values for the dressed molecular tunneling, 
including diagonal tunneling, are provided in Ref. [13] . 

The effective detunings of the dressed molecules are 
evaluated as 



dxdyWt- (x,y)ffW 4 , Q (x,y) 



dK 

VBZ 



(82) 



which is the average of the interacting band structure in 
the BZ. The magnitudes of these detunings, compared 
with the pairing strengths to be discussed, determine 
which of the dressed molecules are nearby the lowest band 
in energy and so must be included dynamically in the 
FRH for a given a s . The Hubbard parameters for the 
pairing are evaluated as 



9i-k,k- 



dxdyWt- (x, y) Hwj (x) w k (x) . (83) 



The open channel fermions couple to the dressed 
molecules only through Feshbach coupling to their closed 
channel components according to the two-channel model, 
see Eqs. ( flifls] ). Thus, 
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(85) 



where we have defined A, 



R, 



R, and used the 



coupling constant at total quasimomentum K, g a j£, is 

l-T/K I rpa 
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In the limit of an infinitely broad resonance where 
oo, we have that 
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(87) 



which depends on the divergent quantity g only para- 
metrically through the derivative of the matrix % K . 
Physically, as the resonance becomes very broad, the 
closed channel population of the eigenstate vanishes, see 
Eq. (77 1. However, the interchannel coupling diverges 



definition Eq. (22 1 in the last line. Hence, the effective 



in this same limit, and the infinite coupling to an in- 
finitesimal closed channel population produces a finite 
limit. In Appendix [A] it is proved that the matrix 
~^X K (-£'k/ Er) /dEjt is positive semidefinite and so the 
pairing is real and the transformation to the FRH has 
the effect of narrowing the resonance irrespective of the 
value of the bare pairing strength g. The behavior of 
the effective pairing Eq. (851 in typical cases is given in 
Ref. [13]. 

We note that, while the FRH is still in general a multi- 
band Hamiltonian with the eigenstate index a as the 
band index, in typical cases we need only a single dressed 
molecular band as compared to the large number of bands 
required for an accurate calculation without dressing the 
closed channel. For example, for the case of a reso- 
nance in the lowest band of the open channel, the lowest 
energy dressed molecule with completely even parity is 
the most relevant to describing the resonance behavior. 
The dressed molecules thus form an ideal dressed closed 
channel, and the FRH is the two-channel lattice model 
between fermions in the lowest band and these dressed 
molecules. Dressed molecules which are off-resonant can 
be adiabatically eliminated from the theory, and their 
effects can be included as a background scattering term 
using e.g. second order perturbation theory. The possibil- 
ity of three- and higher-body effects are not captured by 
the two-channel model, and so such terms do not appear 
in the FRH as written. However, this is a shortcoming of 
the continuum modeling, step 1 from Fig. [I] and can be 
overcome by performing the lattice projection of a more 
complex continuum scattering model. 



VI. CONCLUSIONS 

We have presented a complete two-particle theory 
of two-component fermions in a three-dimensional lat- 
tice potential interacting through a Feshbach resonance 
within the confines of the two-channel model. The two- 
particle lattice bound states are obtained from an eigen- 
value equation which is nonlinear in the energy eigen- 
value. Specializing to the case of a simple cubic lattice, 
we showed that the matrix elements of the interchannel 
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coupling in the Bloch basis have a simple representation 
in terms of the Fourier components of the Bloch basis. 
This Fourier expansion may be found efficiently numer- 
ically by solving a real symmetric eigenvalue problem. 
Analyzing the interchannel overlaps, we demonstrated 
that the effect of an interchannel regularization function 
which has the same symmetries as the lattice may be in- 
terpreted in terms of a Fourier cutoff in a shifted momen- 
tum space which has the symmetries of the momentum- 
space BZ. This result unifies the understanding of the 
renormalization of two-channel theories in the lattice and 
in free space, as renormalization in the latter amounts to 
imposing a spherically symmetric cutoff in momentum 
space. Using the simplified expressions for the interchan- 
nel overlaps, we computed the lattice renormalization at 
zero total quasimomentum, and demonstrated that the 
renormalization may also be interpreted as a cutoff in the 
contributions from higher Bloch bands. We presented a 
scaling relation for the effective closed channel T-matrix 
which allows results for finite regularization cutoff to be 
extrapolated to infinite cutoff. By arranging limits ap- 
propriately, this scaling procedure produces results for 
the bound states which properly accounts for the effects 
of all higher bands. 

The two-particle theory was then applied to a two- 
channel model in which the lowest open channel scat- 
tering band was projected out. We called these bound 
states dressed molecules. When recombined with the low- 
est band of the open channel, the dressed molecules accu- 
rately reproduce the resonance physics at the two-body 
level, including proper renormalization and the effects of 
higher bands. We analyzed the dressed molecular disper- 
sion as a function of s-wave scattering length at zero total 
quasimomentum, and showed that their interpretation as 
an effective closed channel for the two-body resonance is 
valid. We then defined localized Wannier functions from 
these two-body states and gave a phase prescription such 
that the Wannier functions transform with the symme- 
tries of the lattice. Taking matrix elements of the many- 
body Hamiltonian in the basis of the dressed molecu- 
lar Wannier functions, we defined the Fermi resonance 
Hamiltonian (FRH). The FRH has the interpretation of a 
many-body resonance model between unpaired fermions 
in the lowest band and dressed molecules. By symme- 
try and energetic considerations, only a single dressed 
molecule is relevant in typical cases, and so the FRH of- 
ten reduces to an effective two-channel model. The use of 
the full lattice solution to determine the dressed molecu- 
lar Wannier functions not only increases the quantitative 
predictive power of the model versus models which use 
separable approximations to the lattice, but also leads to 
interesting features of the bound states such as tunneling 
along non-principal axes. 

The results in this work can be applied and extended 
to a variety of scenarios. The problem of pairing of 
unequal-mass fermions is relevant to the formation of ul- 
tracold molecules via magneto-association across a Fes- 
hbach resonance [SS]. The interesting case of quasi-low 



dimensional systems with strong s-wave interactions can 
be studied by altering the lattice depths along the princi- 
pal axes in the lattice potential, see Eq. ( 35 1 . In a similar 



vein, our results can be extended to arbitrary lattice ge- 
ometries. In particular, our derivation of the nonlinear 
eigenvalue equation for bound states and our exposition 
of the renormalization in terms of a momentum cutoff 
which has the same shape as the reciprocal-space BZ ap- 
ply to any lattice geometry. Finally, more general inter- 
channel potentials can be accommodated, for example to 
study pairing to higher angular momentum states. 
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Appendix A: Properties of the derivative of the 
effective closed channel T-matrix 



Here we wish to prove that —d^/dE is positive def- 
inite outside of open channel scattering bands. As we 
consider ourselves to be outside of scattering bands, we 
drop the infinitesimal n in Eq. (28). We note that we 
may write —dx^/dE as the Gram matrix 



dE 

by defining the vectors 



/nmq 
sK 



ftK> 



"sK 



(q) 



E 



K 



- E K 

rim 



(q) 



(Al) 



(A2) 



This proves that —dx^/dE is positive semidefinitc. Pos- 
itive definitencss follows from the fact that the set 
{fate} are linearly independent. We prove this by first 
noting that the parameterization Eq. ( 18 ) of functions 
which transform K -periodically as ip (x + a, y + a) = 
e iK ' a V> (x, y) for any Bravais lattice vector a implies that 

( X ) 0nq (x') m ,K-q (y) </>m,K-q (/) 

nmq 

= 5(x-x')5(y-y') (A3) 

within the space of if-periodic functions. That is to say, 
products of Bloch functions with fixed total quasimomen- 
tum form a complete basis for the space of functions 
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which transform with this total quasimomentum under 
translations. This implies that 



lim 



E^"k (qKT* (q) 



(A4) 



nmq 



The elements of h™j£ (q) with s and K fixed, arranged 
as a vector, h S K, thus form an orthogonal set, and so 
the determinant of the matrix with these vectors as rows 
is nonzero. If we now scale the nmq element of each 
vector h S K by the same factor [Ejc — E^ m (q)) , this 
determinant is multiplied by the product of all of these 
factors. Hence, the determinant is nonzero as long as 
each of the factors is nonzero. The factors are nonzero 
so long as the energy 7?k is finite and does not lie in 
one of the scattering bands of the open channel. The 
scaled vectors are the set {f S K:}, and the non-vanishing 
determinant implies that this set is linearly independent, 
completing the proof. 

For the projected model, we instead deal with 



d\ /BE and the relation Eq. (A4) no longer applies, 



as we have excluded the lowest band of the open channel. 
Hence, our results show only that —dj^/dE is positive 
semidefinite. In all numerical cases we have considered, 
the scalar product — T?5 • -Jfer ■ 



OE 



is nonzero. 



Appendix B: Construction and scaling of the 
effective closed channel T-matrix 



The most numerically demanding part of the two- 
particle solution is the construction of the matrix 
X K (E K ) defined in Eq. @. While the overlaps h£g (q) 
may be written as pro ducts of ID functions for the simple 
cubic lattice, see Sec. IV C the construction of x K (-Ek) 
does not factorize due to a nontrivial denominator. In 



this section, we present an overview of efficient numeri- 
cal methods for computing the values of x K (Fk) and for 
performing the scaling Eq. ( 69 ) necessary for a properly 
renormalized theory. 

We begin with a discussion of the integration over the 
BZ in Eq. (28). Because little is known generally about 



the behavior of the integrand for arbitrary K and band 
indices, it is best to look for an adaptive approach which 
refines the integration grid until a desired numerical tol- 
erance is met. For the cubic BZ that we consider, a nat- 
ural choice is the Genz-Malik algorithm [60] , an adaptive 
integration method over hyper-rectangular regions in 3D. 
We define a generator F [J, v] which takes as input a func- 



tion 7 fvl 



and a vector v G 



as the sum of 



function evaluations of I at all unique permutations of 
the elements of v including sign changes. For example, 

F [7, (0,0,0)] = 7 ((0,0,0)) , (Bl) 
F [I, (1, 0, 0)] = / ((1, 0, 0)) + 7 ((0, 1, 0)) + 7 ((0, 0, 1)) 
+ 7 ((-1,0,0)) +7 ((0,-1,0))+/ ((0,0,-1)) . (B2) 



The Genz-Malik algorithm uses the ansatz 
,i ,i ,i 

dx dy dzl (x,y, z) pa Rj (I,w, A) 
-l J-x J-i 

R 7 (I, w,X)= Wl F [I, (0,0,0)} 
+ w 2 F [7, (A 2 , 0, 0)] + w 3 F [I, (As, 0, 0)] 
+ Wi F [7, (A 4 , A 4 , 0)] + w 5 F [I, (A 5 , A 5 , A 5 )] , 



(B3) 



(B4) 



and then chooses the parameters {wi} and {A^} such that 
i?7 is a rule of degree seven, meaning that it exactly inte- 
grates all monomials of degree at most seven |61| . This is 
in analogy to the more familiar Gauss-Legendre quadra- 
ture in ID, where a set of integration points and weights 
are chosen such that the quadrature rule exactly inte- 



grates all polynomials of a given degree 
priate choice of parameters is |60j 

87488 7840 



An appro- 



Wl 



'U'4 



19683 ' 
1600 



w 2 



w 5 



6561 
6859 



w 3 



4960 
19683 



A 



(B5) 
(B6) 



19683 ' ° 19683 ' 

70' A3 = A4= 10' As= 19 

The domain of integration can be deformed to any hyper- 
rectangular region by appropriate scaling of the inte- 
grand. Additionally, as the integration points do not 
include the boundaries, the domain of integration can be 
constructed so as to avoid troublesome function evalu- 
ations for singular integrands. To estimate the error of 
the rule R7 (7, w, A) a rule of degree five, i? 5 (7, w', A), is 
constructed using only function evaluations at the points 
as 

R 5 (I,w',X)=w'- L F[I, (0,0,0)] 



u> 2 F[7,(A 2 ,0,0)] 



^F[7,(A 3 ,0,0)] + <F[7,(A4,A 4 ,0)] 



with 



4456 
243 



980 
243 



140 
729 



w 4 



(B7) 



200 

729 ' 

(B8) 



The algorithm begins by dividing the integration re- 
gion into hyper-rectangular regions, and then using 
Eqs. ( B3 1 and ( |B7[ ) to provide degree 7 and 5 estimates of 
the integral in each region, 7? 7 and i? 5 . The error in each 
region is estimated by the absolute difference |i?7 — R§\, 
and the region with greatest error is identified. We now 
wish to subdivide the region of greatest error into two 
hyper-rectangular regions by dividing the region in half 
along a Cartesian dimension. The dimension to be sub- 
divided is the one with the greatest fourth difference of 
the integrand, estimated as, e.g., 



TV 4) 7(v) 



=(0,0,0) 



(B9) 



=7 ((A 2 , 0, 0)) + 7 ((-A 2 , 0, 0)) - 27 ((0, 0, 0)) 



A 



2 [I ((A 4 , 0, 0)) + 7 ((-A 4 , 0, 0)) - 27 ((0, 0, 0))] . 
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Here, we have placed the center of the region to be 
subdivided at v = (0, 0, 0) for simplicity. We note 
that the fourth derivative estimation uses only points re- 
quired in evaluating Rj and R§. Estimates for the in- 
tegral on the new subregions are provided by applying 
Eqs. (B3| and (B7|, and the process of finding the re- 



gion of largest error, subdividing according to the great- 
est fourth difference, and computing integrals on the new 
subregions is continued until the total error is bounded 
by a user-defined tolerance. A nice feature of the Genz- 
Malik algorithm is that it parallelizes well over subre- 
gions [B3]. Additionally, as discussed in Sec. IV C the 



matrix elements (q) separate along principal axes 
and so the computation of these functions may be cached 
along each Cartesian dimension to avoid overhead in re- 
computation. Similar algorithms can be devised for non- 
rectangular regions [04 to enable the study of the FRH 
on other lattice structures. 

As discussed in Sec. |IVD[ proper computation of 
X K (Ek) first fixes A and converges a summation as the 
number of open channel bands S — > oo before letting 
A — > oo according to the scaling relationship Eq. ( [69|. In 
practice we use shell summation as defined in Sec. IV D| 
to converge x K (Ek) as S — > oo to a desired relative pre- 
cision e re i. The diagonal elements (-Ek) are typically 
larger than the off-diagonal elements, and so it may be 
desirable to put an absolute tolerance e a b s on the off- 
diagonal elements which is comparable to e r ciX^ (Ek) 
but larger than ercixS (Ek) in order to lower the compu- 
tational cost for the same accuracy. That the eigenvalues 
of x K (Ek) are computed to the appropriate relative pre- 
cision with this increased tolerance for the off-diagonal 
elements follows from the fact that the maximal differ- 
ence between eigenvalues of two matrices A and M is 
bounded by the 2-norm difference of the two matrices, 



max |Aj- [A] - Xj [M} \ < \\A - M| 



with Xj [A] the j eigenvalue of A 



(BIO) 



Appendix C: Solution of the nonlinear eigenvalue 
problem 

In this section we focus on the numerical solution of the 



nonlinear eigenvalue equation Eq. ( 27 1 . For simplicity, we 



will focus on the limit of an infinitely broad resonance 
in which g — > oo, hv/E R — > oo, a s — >const such that 
Eq. (27) becomes 



1 



t 



f t (E K /E R )Tl 



k 



(CI) 



This is essentially the limit of a zero-range resonance, 
tb j a — > in which the s-wave scattering length is the 



only relevant parameter. The method here also applies 
to finite width resonances [13]. 



The simplest means to solve Eq. ( CI ) is to fix the eigen- 



value E'k and then solve Eq. (CI) as a linear eigenvalue 
problem for a s . We will refer to this method of finding 
solutions as the linear method. The computation of Hub- 
bard parameters requires the bound state solution in the 
entire total quasimomentum BZ at fixed a s rather than 
fixed Ek, see Eq. (851, for example. 

To solve Eq. (CI I for fixed a s , we first restate the prob- 
lem as as 



X* (E K /E R ) 



7T 1 



-K 



(C2) 



We assume that we have some guess at the eigenen- 
ergy i?K, which will not necessarily be the solution but 
will sit on some continuous path traced by an eigenvalue 
Xica s u (Ek) of T K;as (Ek) provided that the guess does 
not lie in one of the scattering bands of the open chan- 
nel [66] . Hence, we can perform a root-finding method on 
AKn,a (Ek) to search for exact eigentuples. The Newton 
iteration yields the shift in the approximate eigenvalue 
Ek/Er at the i th iteration as 



T?-T, 



k 



a s (Ek) 



T 



k 



YK. T 



'K;5 3 (Ek) 



(C3) 



where Y^- is the eigenvector of Tk ; 5 s (EjcJ which is 
smoothly connected to our initial guess. As a s is a fixed 
parameter appearing in Tks s (Ek), the denominator is 
in fact the quantity • \' (E^_/E R ) ■ appearing 
in the normalization factor of the two-particle solution. 
One updates the guess at the eigenenergy using this shift 
and iteration continues until the shift drops below a given 
tolerance. This Newton-Raphson method can be proven 
to be quadratically convergent locally [66] . 

The complete procedure for finding the band structure 
at fixed a s is first to find the exact solutions (a s , Ek, Y k ) 
at the extremal points of the BZ where all elements of 
K are either or —it /a. This is done by choosing a 
range of -Ek and solving the associated linear eigenprob- 
lem for d s by the linear method. As discussed in Sec. |VB[ 
these solutions can be classified according to their par- 
ity unambiguously. Now, a s is fixed and the solutions at 
the extremal points are used as guesses for solutions at 
fixed a s and K with a desired parity. A nice feature of 



the Newton-Raphson iteration Eq. (C3) is that we can 
search along directions which are smoothly connected to 
guesses with definite parity, and so the Kohn construc- 
tion discussed in Sec. IV Bl can be enforced. 
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